;+
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;
; Read MM5 Output file and extract data
;
; Usage: x=read_mm5(file, varlist, /onefield)
;
; By claydodo
; Based on the matlab script version of Klaus Wyser
;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;-

function read_mm5, filein, varlist, onefield=onefield, maxstep=maxstep

	if (n_elements(varlist) eq 0) then begin
	;;; .m: varlist={};
	endif else begin
		varlist=strlowcase(varlist)
	endelse

	if (n_elements(onefield) eq 0) then onefield=0
	
	if (n_elements(maxstep) eq 0) then maxstep=25

	if (n_elements(filein) eq 0) then filein='fort.10'

	nvar=0L ;;; .m: varfound={}, but idl cannot make zero element [].
	istep=0L

	openr, fid, filein, /get_lun, /swap_if_little_endian
	flag=0L
	
	dum=0UL ;;; should use F77_UNFORMATTED in openr instead.
	bhi=lonarr(50,20) ;; 
	bhr=lonarr(20,20)
	bhic=bytarr(80,20,50) ;;
	bhrc=bytarr(80,20,20) ;;
	ndim=0L ;long
	start_index=lonarr(4)
	end_index=lonarr(4)
	xtime=0.0 ;float
	staggering=bytarr(4)
	ordering=bytarr(4)
	current_date=bytarr(24)
	name=bytarr(9)
	unit=bytarr(25)
	description=bytarr(46)



	while ~ eof(fid) do begin
		readu, fid, dum		
		readu, fid, flag
		readu, fid, dum

		case flag of
		0L : begin
			print, '...found big header'

			;read big header
			readu, fid, dum
			readu, fid, bhi
			readu, fid, bhr
			readu, fid, bhic ;;.m did something more here
			readu, fid, bhrc
			readu, fid, dum

			if (n_elements(varlist) eq 0) then begin
				print, 'File type BHI(1,1)=', string(bhi(0,0)) ;;!start with 0!
			endif

		end
		1L : begin
			;read subheader
			readu, fid, dum
			readu, fid, ndim
			readu, fid, start_index
			readu, fid, end_index
			readu, fid, xtime
			readu, fid, staggering
			readu, fid, ordering
			readu, fid, current_date
			readu, fid, name
			readu, fid, unit
			readu, fid, description
			readu, fid, dum

			; see if variable exists
			varname=strtrim(strlowcase(name));
			varname=strjoin(strsplit(varname,/extract),'_')
			;print, varname
			if (n_elements(varfound) eq 0) then begin
				condition2 = 1 
			endif else begin
				condition2 = ( total(strcmp(varname, varfound)) eq 0 )
			endelse

			if (total(strcmp(varname, varlist, /fold_case)) gt 0 &&   $
				condition2 ) then begin ;;; varfound!

				case ndim of
				1 : begin
					exeres=execute(varname+'=fltarr(end_index[0]-start_index[0]+1, maxstep)')
				end
				2 : begin
					if (strcmp(strtrim(ordering),'CA')) then begin
						exeres=execute(varname+'=fltarr(end_index[0]-start_index[0]+1, end_index[1]-start_index[1]+1)')
					endif else begin
						exeres=execute(varname+'=fltarr(end_index[0]-start_index[0]+1, end_index[1]-start_index[1]+1, maxstep)')
					endelse
				end
				3 : begin
					exeres=execute(varname+'=fltarr(end_index[0]-start_index[0]+1, end_index[1]-start_index[1]+1, end_index[2]-start_index[2]+1, maxstep)')
				end
				else : begin
					print, 'unknown ndim=', string(ndim)
					stop
				end
				endcase

				nvar++
				if (nvar eq 1 ) then begin
					varfound=varname
					varorder=strtrim(ordering)
					varstag=strtrim(staggering)
				endif else begin
					varfound=[varfound,varname]
					varorder=[varorder,strtrim(ordering)]
					varstag=[varstag,strtrim(staggering)]
				endelse
			endif else begin
				if (n_elements(varlist) eq 0) then begin
					print, strtrim(name), strtrim(unit), strtrim(description), end_index-start_index+1
					;;; .m has format here
				endif
			endelse

			datadim=end_index-start_index+1
			;print, datadim
			readu, fid, dum
			data=fltarr(long(product(datadim)))
			readu, fid, data
			readu, fid, dum
			if (total(strcmp(varname, varlist, /fold_case)) gt 0 ) then begin
				print, '......found ', varname, ' at time ', string(xtime)

				case ndim of
				1 : begin
					data=reform(data,datadim[0],1)
					exeres=execute(varname+'[*,istep]=data')
				end
				2 : begin
					if (strcmp(strtrim(ordering), 'CA')) then begin
						data=reform(data, datadim[0], datadim[1])
						exeres=execute(varname+'[*,*]=data')
					endif else begin
						data=reform(data, datadim[0], datadim[1])
						exeres=execute(varname+'[*,*,istep]=data')
					endelse
				end
				3 : begin
					data=reform(data, datadim[0],datadim[1], datadim[2])
					exeres=execute(varname+'[*,*,*,istep]=data')
				end
				else : begin
					print, 'unknown ndim=', string(ndim)
					print, '*** RETURN 2 ***'
					stop
				end
				endcase

			endif
		end
		2L : begin
			if (istep eq 0) then begin
				ftime=xtime
			endif else begin
				ftime=[ftime,xtime]
			endelse
			istep++
			print, '...end of time period ', string(xtime)
			if (n_elements(varlist) eq 0) then return, 0
			if (onefield eq 1) then break
		end
		else : begin
			print, 'no valid flag: ', string(flag)
		end
		endcase
		
	endwhile
	
	print, 'end of file'
	free_lun, fid
	; remove unused part of variables
	istep=max([0, istep-1])
	for i=0, nvar-1 do begin
		exeres=execute('x='+varfound[i]);
		siz=size(x)
		;print, siz
		if ( varstag[i] eq 'C' ) then begin ;staggering
			siz[1] -= 1
			siz[2] -= 1
		endif
		case size(x, /n_dimensions) of
		2 : begin
			if ( ~ strcmp(varorder[i],'CA') ) then begin
				exeres=execute(varfound[i]+'=x[0:siz[1]-1, 0:istep]')
			end
		end
		3 : begin
			exeres=execute(varfound[i]+'=x[0:siz[1]-1, 0:siz[2]-1, 0:istep]')
		end
		4 : begin
			exeres=execute(varfound[i]+'=x[0:siz[1]-1, 0:siz[2]-1, *, 0:istep]')
		end
		else : begin
			print, '*** RETURN 3 ***'
		end
		endcase
	endfor

	; start filling output structure
	totalcommand='x={xtime:ftime'
	for i=0, n_elements(varlist)-1 do begin
		if (total(strcmp(varlist[i], varfound)) eq 0) then begin
			print, '*** Warning: ', varlist[i], ' not found ***'
		endif else begin
			totalcommand = totalcommand +','+varlist[i]+':'+varlist[i]
			;exeres=execute('x={x,'+varlist[i]+':'+varlist[i]+'}')
		endelse
	endfor
	totalcommand += '}'
	exeres=execute(totalcommand)

	return, x

end
